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ABSTRACT 

We identify a new hydrodynamical instability in protoplanetary discs that may arise due to 
variations in the dust-to-gas ratio and may lead to concentration of dust grains within a disc. 
The instability can arise due to dust settling, which produces a vertical compositional entropy 
gradient. The entropy gradient drives a baroclinic instability that is capable of creating toroidal 
gas vortices that gather dust into rings. Such dust rings are potentially observable via contin¬ 
uum emission of the dust or scattered light. Indeed, this instability may offer an explanation 
for the rings recently observed in the discs around the young stars HL Tau and TW Hya that 
does not rely on clearing by protoplanets. The instability may also have wider ramifications, 
potentially aiding dust agglomeration, altering the radial migration of larger planetesimals, 
and modifying angular momentum transport within a disc. 
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1 INTRODUCTION 

During the early phases of planet formation, dust grains are hy¬ 
pothesised to agglomerate into planetesimals. However, collisions 
between planetesimals larger than a few centimetres in size are gen¬ 
erally thought to lead instead to bouncing or frag mentation (|Wei- 
denschilling|1977{|Weidenschilling & Cuzzi|1993||Blum & Wurm 

2000). In addition, in laminar gas discs, metre-sized planetesimals 
are predicted to suffer rapid orbital decay due to gas drag ( |Wei-| 
denschilling 1977). Various mechanisms have been proposed to 
concentrate dust and planetesimals at particular locations in discs, 
either to avoid radial migration or to allow them to grow (e.g. 


Safronov 

1969, 

Goldreich & Ward 1973, Nakagawa et al. 

1986 

Weidensc 

lilling 

1980| Barge & Sommeria 1995 ; Tanga et al. 

1996 

Rice et a 

200‘ 

, jYoudin & Goodman 2005, Johansen & Klahr 


2005| ). In particular, baroclinc instabilities in protoplanetary discs 
may form large-scale vortices (j Klahr & Bodenheimer|2003[ | Klahr | 
2004, Petersen et al. 2007), and these have been suggested as a pos- 


sible mechanism to trap dust particles jBarge & Sommeria 

1995, 

Tanga et al. 1996, Godon & Livio 1999; Johansen et ak]2004 

Klahr 

& Bodenheimer|2006). In these past studies, the baroclinic insta- 


bilities have been proposed to arise due to radial thermal entropy 
gradients in the disc and produce vortices coplanar with the disc. 


A dust particle moving through gas experiences drag. The 
form of the drag force may vary considerably as a function of the 
grain and gas properties fWeidenschilling| 1977] ), but if the mean 
free path of the gas molecules is larger than the dust particle’s 
radius, the characteristic timescale for the dust particle’s speed 
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relative to the gas to decay (the stopping time) depends linearly on 
its density and size, and inversely on the density of the gas. Thus, 
the behaviour of dust particles in a gaseous disc around a young 
star depends on their size. Small particles (<C 1 mm for typical 
protoplanetary disc densities) have short stopping times and are 
well-coupled to the gas. The dust-gas mixture behaves almost 
as a single fluid, with a lower pressure for a given density than 
for pure gas (since the dust does not contribute to the pressure). 
Large particles (typically 10 m) are poorly coupled and move 
on almost ballistic trajectories, barely affecting the gas. Particles 
with intermediate sizes affect the inertia of the dust-gas fluid, but 
also migrate relative to the gas. In particular, if such particles 
are initially well-mixed in the gaseous disc surrounding a young 
star, they will tend to slowly settle (Gold reich & Ward||1973| > 
towards the mid-plane of the disc under the action of the vertical 
component of the gravitational field exerted by the star. 

Dust settling is a phase separation process. One therefore ex¬ 
pects to have a strong entropy gradient at the interface between 
the settled, gas/dust mixture, and the purer gas phase above. In a 
gaseous protoplanetary disc surrounding a young star, settling of 
dust towards the mid-plane of the disc will produce a vertical gradi¬ 
ent of the dust-to-gas ratio above the disc’s mid-plane. In this letter, 
we show that this vertical compositional entropy gradient may pro¬ 
duce an azimuthal baroclinic instability that drives toroidal vortices 
within the disc. In Sectionp] we investigate analytically the criteria 
for instability. In SectionR] we describe the numerical method and 
initial conditions that we use for hydrodynamical simulations of 
dusty protoplanetary discs. Results from our calculations are pre¬ 
sented in Section]?] and in Section[5]we give our conclusions. 
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Figure 1. The time evolution of the azimuthally-averaged dust-to-gas ratio 
(colour rendering, maximum e ~ 0.025) and gas density (greyscale, rang¬ 
ing from log 10 (pgas) = [—16, —11]) for a locally-isothermal gaseous cir- 
cumstellar disc with H/r = 0.05 and an initial dust-to-gas ratio, e = 0.01 
comprised of 1 mm dust grains. The dust initially settles toward the mid¬ 
plane (top panel). After the dust reaches its quasi-stationary state, the com¬ 
positional baroclinic instability rapidly develops, first at small radii (middle 
panel), and then throughout the disc (lower panel). The radial wavelength of 
the instability is approximately equal to the thickness of the quasi-stationary 
dust layer. The calculation was performed with N p = 3 x 10 6 particles. 


2 HYDRODYNAMICAL INSTABILITY 


Baroclinity is described mathematically by the curl of the pressure 
gradient term in the fluid momentum equation 

V x ( _ ~vpj =-T V/ o X VP, (1) 

where P and p are the fluid pressure and density, respectively. In 
an axisymmetric dusty gas in which the dust is well coupled to the 
gas and the mixture behaves as a single fluid, the only non-zero 
component of equation [T] in cylindrical coordinates is given by 


1 P \dp dK K (dp de 

p 1 p Kp 2 [dz dr 1 l + e\dr dz 


dp de 
dz dr 


5 

( 2 ) 


where ee is the azimuthal unit vector and, for generality, we have 
assumed a barotropic equation of state given by P = KpJ as = 

K , where K(r), the total density of the fluid p = p gas + 

Pdust is the sum of the densities of the gas and dust, and e = 
Pdust/pgas is the dust-to-gas ratio. A non-zero baroclinity will 
change the vorticity, potentially producing vortices. However, the 
presence of baroclinity by itself does not necessarily lead to insta¬ 


bility due to the stabilising effect of rotation. In a typical gaseous 
accretion disc, the vertical density stratification (first term on the 
right hand side of equation^ is not strong enough to make the disc 
unstable. However, as can be seen in the last term on the right hand 
side of equation [ 2 ] the presence of a compositional gradient may 
significantly increase baroclinity, potentially making a disc unsta¬ 
ble. As has been beautifully explained in the context of thermally 
driven baroclinic instabilities (e.g. |Lyra & Klahr||2011| ), by anal¬ 
ogy, a compositional instability may be initiated due to the presence 
of strong gradients in the dust-to-gas ratio, and the resulting baro¬ 
clinity will generate vorticity. The emerging vortices mix the fluid, 
moving high entropy well-mixed fluid from the dusty settled layers 
into the lower entropy purer gas, thus reinforcing the instability. 

The two Solberg-Hpiland stability conditions for a gaseous 
disc ( |Tassoul|2000| ) are 

4^- -VP- VS> 0, 

r 3 dr Cpp ( 3 ) 

dS/dz > 0, 


where j = r 2 9 is the Keplerian specific angular momentum and Q 
is the Keplerian angular frequency, c p is the specific heat capacity 
of the gas at constant pressure, and S is the entropy of the fluid. 
These equations will still be valid in a dusty gas as long as the mix¬ 
ture behaves as a single fluid (i.e. the dust stopping time is short) 
and the specific heat capacity of the mixture is dominated by that 
of the gas (i.e. the dust-to-gas ratio is small). 

The first of the Solberg-Hpiland conditions will be violated if 
the Brunt-Vaisala frequency squared (the second term) dominates 
over the epicyclic frequency squared (the first term). Due to dust 
settling, the vertical gradient in the dust-to-gas ratio will initially 
dominate the radial gradient, so \VS\ ~ \dS/dz\ and 
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where ( |Cabbot| 1984] ), 
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d\ogp 
d log P 


1 
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(5) 


c v is the specific heat capacity of the gas at constant volume, the 
sound speed is c s = y/jP/p, and g z — — fl 2 z is the vertical com¬ 
ponent of the acceleration due to the star’s gravity. The quantity in 
equation [5] is usually zero in an equilibrium situation, but in a con- 
vectively stable region typically adopts values AV ~ 0.1 ( |Cab-| 
|bot|l984| . In our locally-isothermal disc models discussed below, 
AV = 0, so the vertical entropy gradient arises purely due to the 
dust-to-gas ratio vertical gradient. We can then write 
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Therefore, the system is potentially unstable at 2 > 0 if 
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Figure 2. For the same case as shown in the lower panel of Fig. 1, we plot the velocity field (vectors) and the azimuthally-averaged dust density (colour 
rendering) and gas density (greyscale, ranging from log 10 (p ga s) = [—17, —12] in g cm -3 ). The toroidal vortices within the dust layer are clearly evident, 
and the dust density is highest at the centres of the vortices (i.e. producing concentric rings). 


where we have defined the vertical scale height of the disc H = 
c s /Q. The former equation implies that it becomes harder for 
the instability criterion to be satisfied for small z, and the most 
favourable location for the development of the instability is 

z _ 1 _^ de 

H ~ VtAV ^ dz 

If AV = 0, the criterion for instability is more easily satisfied for 
increasing z. 


2(l + e) VtAV 
H 


( 8 ) 


3 NUMERICAL METHOD 

To perform numerical simulations of dusty protoplanetary discs, we 
have used the three-dimensional SPH ( |Lucy|1977||Gingold & Mon-] 
| aghan ||T977] ) code, sphNG (Ayliffe et al.||201l| ), which has been 
used to study a variety of astrophysical fluid dynamical systems. 
For the simulations performed for this letter, the code uses two 
populations of SPH particle: gas and dust. The latter only evolves 
under the action of gravity and a drag force. In order to perform the 
time-integration of the dust-gas drag force, use the semi-implicit 
integration method of |Loren-Aguilar & Bate|p0l4) . 

The simulations performed followed the time evolution of 
locally-isothermal (T ~ r -1 and independent of time) circumstel- 
lar gas/dust discs. The central star was modelled as a gravitational 
point mass, M* = 1 M©, the mass of the gas component between 
the inner and outer radii was set to be « 0.003 M© (similar to the 
minimum-mass Solar nebula), and the mass of the dust component 
was determined by the initial dust-to-gas ratio. The gaseous com¬ 
ponent of the disc had an outer radius of r = 16.1 AU, with an 
inner boundary at r = 2.6 AU, chosen to avoid having particles 
with excessively short integration time-steps. We modelled discs 
with vertical scale heights of H/r = 0.025, 0.05, or 0.10 and sur¬ 
face density profiles E oc r _1 ^ 2 . Self-gravity was neglected. We 
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performed calculations with dust-to-gas ratios of e = 0.001, 0.01, 
or 0.1. Calculations were typically performed with resolutions of 
iV p = 6 x 10 5 SPH particles, but some cases were performed with 
N p = 3 x 10 6 . In each case, approximately 1/3 of the particles were 
dust particles. The initial dust-to-gas ratios were small enough to 
avoid significant perturbations of the gas hydrostatic equilibrium. 
After gas component relaxation, dust grains were introduced into 
the disc, following the gas density distribution, but with minimum 
and maximum orbital radii of r — 5.2 AU and r m 15.6 AU, 
respectively. We performed simulations assuming spherical dust 
grains with radii of 1 gm, 0.1 mm, 1 mm, 1 cm, 1 m and 10 m, 
each with intrinsic dust grain densities of 3.0 g cm -3 . 

Calculations performed at different numerical resolutions pro¬ 
duced qualitatively similar results in terms of whether the insta¬ 
bility developed or not, though the radial spacing of the vortices 
tended to decrease somewhat as the resolution was increased. We 
performed calculations using different SPH kernels: the M4 cubic, 
the M6 quinti c (|Monaghan & Lattanzio|1985| ), and the C4 Wend- 
land kernel jWcndland| 1 995f ~ Very similar results were obtained 
with different kernels. The figures in this letter were all produced 
from a calculation using N p = 3 x 10 6 and the M4 kernel. 


4 RESULTS AND DISCUSSION 

To study the instability, we performed a series of three-dimensional 
smoothed particle hydrodynamic (SPH) simulations of dust settling 
in gaseous locally-isothermal Keplerian discs. As a typical case, 
consider a gas disc with H/r — 0.05 with a grain size of 1 mm 
and a dust-to-gas ratio of 0.01. For radii of 5-10 AU, the dust set¬ 
tling times are of order a hundred years for such partially coupled 
grains. The dust particles reach a stationary solution where their 
vertical velocity becomes negligible jNakagawa et al.|1986) , pro¬ 
ducing a dust layer centred on the mid-plane with a width inversely 
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Figure 3. Surface density renderings of the gas (top) and dust (bottom) after 
the instability has developed throughout the disc, for the same calculation 
as that depicted in Figs. 1 and 2. The same calculation as that shown in 
Figs. 1 and 2, but after the instability has developed throughout the disc. 
The dust is conglomerated into concentric rings with peak densities located 
at the centres of the toroidal vortices within the dust layer. 



Figure 4. A volume density rendering of the dust component whose surface 
density is shown in Fig. 3. 


proportional to their stopping time. As can be seen in Fig. |T| soon 
after reaching the quasi-stationary phase, the disc starts to develop 
the baroclinic instability. Toroidal vortices develop within the dust 
layer (Fig. [2), creating a sustained convective motion (see the video 
at http://www.astro.ex.ac.uk/people/mbate/Animations/dustl.html 
for an illustration of the instability development). Dust particles 
collect at the centres of the vortices, so that looking down on the 
disc the dust conglomerates into rings (lower panel of Fig. [3j. 

We find that the radial spacing of the vortices (and therefore 
the dust rings) is approximately equal to the thickness of the dust 
layer when the instability begins. Therefore, the radial spacing is 
expected to increase for thicker and/or more massive gaseous discs, 


or smaller or less dense dust grains. If the radial spacing between 
rings is large enough, they may be observable at sub-millimetre 
and millimetre wavelengths. Furthermore, in Figs. [2] and |4] it can 
be seen that the dust layer is corrugated, with the vertical peaks 
being out of phase with the density maxima of the dust rings. Full 
radiative transfer modelling is beyond the scope of this letter, but 
we note that because of the vertical corrugation of the dust layer the 
rings may also be visible in scattered light at optical wavelengths or 
due to thermal shadowing at near-infrared wavelengths. Indeed, it 
is tempting to speculate that the rings seen in the scattered light im¬ 
ages of the TW Hya disc ( |Debes et al.|2013| ) or the recent Atacama 
Large Millimetre Array (ALMA) observations of the HL Tau disc 
which displays several concentric rings ( [Brogan et al.|2015) may be 
manifestations of this instability. Since self-gravity is neglected, the 
calculations are scale-free as long as dimensionless quantities are 
kept constant. For example, scaling the simulations to cover radii 
eight times larger (outer radius 128 AU - the approximate radius 
of the HL Tau disc), to maintain the same dimensionless stopping 
time, either the gas mass of the disc must be increased by a factor 
of 64, or the size of the dust particles must be reduced the same 
factor (or a combination of the two). Other mechanisms that may 
produce dust rings in discs include embedded protoplanets (e.g. 

| Paar dekooper & Mellema 2 004]), o r the photoelectric instability 
proposed by [Lyra & Kuchner| ( |2013| ) in low density discs. 

From our other calculations with different parameters, we find 
that the higher the dust-to-gas ratio, the more easily and quickly 
the instability develops. However, we emphasise that the instabil¬ 
ity does not require large dust-to-gas ratios. The exact behaviour 
depends on the stopping time, which in turn depends on the dust 
grain size and density, and gas density. We find that grains with 
sizes smaller than 0.1 mm are tightly coupled to the gas, and the 
disc does not achieve a sufficient compositional gradient to drive 
the instability, at least over the timescales we have modelled. Con¬ 
versely, grains with sizes of > 10 cm or more are too weakly cou¬ 
pled to the gas for it to act as a single fluid and drive the instability. 
They rapidly settle into very thin layers near the disc mid-plane 
without any evidence of the instability. The exact cut off depends 
on the dust-to-gas ratio and the gas density. For our low-density 
discs with gas scale-heights of H/r = 0.10, grains with sizes of 
1 cm are too weakly coupled to the gas to drive the instability even 
with an initial dust-to-gas ratio of 0.1. With H/r = 0.05, the in¬ 
stability develops with 1 cm grains and a dust-to-gas ratio of 0.1, 
but not 0.01. With H/r — 0.025, 1 cm grains drive instability for 
dust-to-gas ratios of 0.01 and higher. 

These numerical results are in good agreement with the an¬ 
alytical model described above, since the closer to the mid-plane 
the dust is (i.e. the larger the grains, or lower the gas density), the 
harder it should be for the instability to develop. In equation [6] one 
can see that the contribution of the vertical gradient of the dust- 
to-gas ratio to the Brunt-Vaisala frequency is proportional to the 
product of the entropy gradient induced by it, and g z . Hence, for z 
values close to the mid-plane, a very strong gradient of the dust-to- 
gas ratio will be necessary to compensate for the weakness of the 
gravitational field. Close to the mid-plane, the instability will not 
develop unless \de/dz\ > 1 /z (equation [7}. 

In past studies, baroclinic instabilities in discs have been pro¬ 
posed to arise due to radial thermal entropy gradients ( [Klahr &| 
|Bodenheimer||2003| |Klahr||2004| |Petersen et al.||2007) . Not only 
is the driving mechanism completely different from the one we dis¬ 
cuss in this letter, but radial entropy gradients produce vortices that 
are coplanar with the disc rather than toroidal. Some recent local 
( |Bai & Stone|2010| ) and global ( |Kowalik et al.|2013) simulations of 


© ??? RAS, MNRAS 000, 
























Toroidal vortices in protoplanetary discs 5 


dusty discs have investigated the streaming instability fYoudin &| 
|Goodman|2005| | Johansen & Youdin|2007a|b| > and found ring-like 
features. However, |Kowalik et al.| ( 2013) exclude the vertical grav¬ 
itational acceleration so there is no dust settling, they do not obtain 
complete rings, and they find that the perturbation wavelength de¬ 
creases with increasing dust-to-gas ratio whereas we find no such 
dependence. |Bai & Stone[ ( [2010) included vertical gravity and stud¬ 
ied the stability of very thin dust layers to Kelvin-Helmholtz and 
streaming instabilities. Although they performed local simulations 
with much thinner dust layers than ours (of order 1/100 the gas 
scale height), in some of their calculations they find radial os¬ 
cillations of the dust density similar to what we find here. They 
ruled out Kelvin-Helmholtz instability as the driving mechanism 
due to the absence of structure in the azimuthal direction (as we 
do), and attributed the radial structures as ‘most likely due to’ the 
streaming instability. Instead, we propose that compositional en¬ 
tropy gradients due to the varying dust-to-gas ratio are crucial for 
the formation of these axisymmetric rings. The instability in Fig. 
[l] takes only « 12 orbital periods to develop, whereas the growth 
time of the streaming instability is expected to be « 60 — 200 
orbital periods |Youdin & Goodman§2005} for these parameters 
(Pdust/Pgas = 0.02, dimensionless stopping times rs ~ 0.02 at 
the mid-plane and rs ~ 0.05 at the edge of the dust layer), which 
we note is comparable to the dust settling time. Further investiga¬ 
tion is certainly warranted. 


5 CONCLUSIONS 

We have shown that dusty protoplanetary discs may be unstable 
to a baroclinic instability that arises due to vertical compositional 
entropy gradients that are naturally produced during dust settling. 
The baroclinic instability produces toroidal gas vortices that gather 
the dust into rings. These rings are potentially observable at (sub- 
millimetre wavelengths using ALMA, and with high-resolution 
telescopes operating at optical or near-infrared wavelengths. The 
susceptibility of dusty protoplanetary discs to develop toroidal 
vortices has implications for a wide range of disc processes that 
warrant future investigation, but are beyond the scope of this letter. 
It is unclear how the compositional baroclinic instability may 
interact with other instabilities in protoplanetary discs, such as 
the dust streaming instability, the thermal baroclinic instability, 
the magneto-rotational instability (Bal bus & Haw ley 199 lj, or 
gravitational instability. In the calculations presented in this letter, 
due to the axisymmetric initial conditions, the baroclinic instability 
produces global toroidal vortices. Although compositional entropy 
gradients will be expected to develop in more structured discs, the 
manifestation of the resulting baroclinic instabilites may be differ¬ 
ent, and they may also alter angular momentum transport within 
the disc. Grains larger than ~ 10 cm cannot drive the instability 
because they are insufficiently coupled to the gas, however, their 
radial migration may be significantly altered by the gas vortices 
produced due to small grains. If the migration of metre-sized 
objects can be slowed or stopped this may aid planet formation. 
Finally, the long-term evolution of the vortices and dust rings 
needs to be examined, and whether or not the conglomeration of 
dust within the toroidal vortices may promote planetesimal growth, 
either by collisional agglomeration or gravitational instability. 

ACKNOWLEDGMENTS 

The figures were created using SPLASH ( [Price | 
|2007| >, a SPH visualization tool publicly available at 
http://users.monash.edu.au/^dprice/splash 


This work was supported by the STFC consolidated grant 
ST/J001627/1, and by the European Research Council un¬ 
der the European Community’s Seventh Framework Programme 
(FP7/2007-2013 grant agreement no. 339248). This work used the 
DiRAC Complexity system, operated by the University of Leices¬ 
ter IT Services, which forms part of the STFC DiRAC HPC Facility 
(www.dirac.ac.uk). This equipment is funded by BIS National E- 
Infrastructure capital grant ST/K000373/1 and STFC DiRAC Op¬ 
erations grant ST/K0003259/1. DiRAC is part of the National E- 
Infrastructure. This work also used the University of Exeter Su¬ 
percomputer, a DiRAC Facility jointly funded by STFC, the Large 
Facilities Capital Fund of BIS and the University of Exeter. 

REFERENCES 

Ayliffe B. A., Laibe G., Price D. J., Bate M. R., 2011, MNRAS, 
423,1450 

Bai, X.-N., Stone J.M., 2012, ApJ, 722, 1437 

Balbus S. A., Hawley J. F., 1991, ApJL, 376, L214 

Barge P., Sommeria J., 1995, A&A, 295, LI 

Blum, J., Wurm, G., 2000, Icarus, 143, 138 

Brogan C. L. et al. (ALMA Partnership), 2015, ApJL in press 

Cabbot W., 1984, ApJ, 277, 806 

Debes J. H., Jang-Condell H., Weinberger A. J., Roberge A., 
Schneider G., 2013, ApJ, 771, 45 
Gingold R.A., Monaghan J.J., 1977, MNRAS, 181, 375 
Godon R, Livio M., 1999, ApJ, 523, 350 
Goldreich P., Ward W. R., 1973, ApJ, 183, 1051 
Johansen A., Andersen A.C., Brandenburg A., 2004, A&A, 417, 
361 

Johansen A., Klahr H., 2005, ApJ, 634, 1353 
Johansen A., Youdin A., 2007, ApJ, 662, 613 
Johansen A., Youdin A., 2007, ApJ, 662, 627 
Klahr H. H., 2004, ApJ, 606, 1070 
Klahr H. H., Bodenheimer P., 2003, ApJ, 582, 869 
Klahr H., Bodenheimer P., 2006, ApJ, 639, 432 
Kowalik K., Hanasz M., Woltanski D., Gawryszczak A., 2013, 
MNRAS, 434, 1460 

Loren-Aguilar P., Bate M. R., 2014, MNRAS, 443, 927 
Lucy L. B., 1977, AJ, 82, 1013 
Lyra W., Klahr H., 2011, A&A, 527, 138 
Lyra W., Kuchner M.,, 2013, Nature, 499, 184 
Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135 
Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375 
Paardekooper S.-J., Mellema G., 2004, A&A, 425, L9 
Petersen M. R., Julien K., Stewart G. R., 2007, ApJ, 658, 1236 
Price D., 2007, Publ. Astron. Soc. Australia, 24, 159 
Rice W. K. M., Lodato G., Pringle J. .E., Armitage P. J., Bonnell 
I. A., 2004, MNRAS, 355, 543 

Safronov Y. S., 1969, Trans. NASA TTF 677, 1972. Nauka Press, 
Moscow 

Tanga R, Babiano A., Dubrulle B., Provenzale A., 1996, Icarus, 
121, 158 

Tassoul J.-L., 2000, Stellar Rotation. Cambridge Univ. Press, 
Cambridge 

Weidenschilling S. J., 1977, MNRAS, 180, 57 
Weidenschilling S. J., 1980, Icarus, 44, 172 
Weidenschilling S. J., Cuzzi J. N., 1993, Protostars and Planets 
III. University of Arizona Press, Tucson, 1031 
Wendland H., 1995, Adv. Comp. Math., 4, 389 
Youdin A. .N, Goodman J., 2005, ApJ, 620, 459 


© ??? RAS, MNRAS 000, [IJ^ 


















